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Abstract 

Flux pinning is still the main limiting factor for the critical current of the high-temperature 
superconductors in high fields. In this paper we model the field dependence of the critical current 
in thin films with columnar defects aligned with the field. The characteristic shape of the critical 
current of the BaZrC>3 doped YBa2Cu30"6+x thin films is reproduced and explained. The model 
is based on solving the Ginzburg- Landau equations with columnar defects present in the lattice. 
The size of the columnar defect is found to be of key importance in explaining the rounded shape 
of the critical current of the BaZrC>3 doped YBa2Cu306+x thin films. It is also found that the size 
of the rod changes the long range order of the vortex lattice. 

PACS numbers: 74.25.Ha, 74.25.Wx, 74.62.Dh 
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I. INTRODUCTION 



The critical current density, j c , of high-temperature superconductors in magnetic field is 
mainly limited by flux pinning^. Most of the envisioned applications, such as generators and 
fault current limiters, require the superconductors to have high critical current at magnetic 
fields of 3-5 T. In this range j c is proportional to the fraction of vortices pinned in different 
pinning sites^. Each individual pinning site provides a pinning force, f p , which depends on 
the type of the pinning site. Typical pinning sites in YBa 2 Cu306+ x (YBCO) superconductors 
are dislocations, twins, antiphase boundaries, impurities and grain boundaries^. Out of these 
dislocations have been found most effective^, since they have the same shape as the vortices. 
Their effectiveness is limited by the small size of the core of the dislocation, 0.3 nm^, which 
is only a fraction of the vortex size in YBCO (£ = 1.5 nm at K-). Twin planes pin 
vortices effectively if the current is parallel to the twin plane and thus the Lorentz force 
perpendicular to the plane, unfortunately if the current is perpendicular to the twin plane, 
the plane channels the vortices for easy movement, and therefore the j c decreases^. 

Doping the YBCO with e.g. BaZr03 (BZO) or BaSnOa (BSO) leads to formation of 
nonsuperconducting nanorods into the superconducting matrix- - — . The nanorods act as 
very efficient pinning sites, since their diameter is around 5-10 nmii^ and therefore around 
the same size as the vortex core. The nanorods are c-axis oriented and increase j c at high 
fields especially when B || q^^. 

The shape of the j c {B) is most often described with the accommodation field, B*, where 
the low field plateau ends and the exponent a, which describes the decrease of j c above 
B* with field £- Q &±±^ri&. For typical undoped YBCO thin films the B*= 40-100 mT and 
a ~ 0.6^. In BZO-doped samples the B* increases up to 0.5 T and a decreases to 0.2 
- 0.4M-. The a-value observed in undoped films is predicted by theories of strong sparse 
pinning sites^- 1 ^, but the lower value in the doped films has not been predicted, and finding 
a simple explanation for the lowered a value seems to be difficult due to the vortex-vortex 
interactions involved. 

The problem in modelling flux pinning is that all the real samples contain many different 
kinds of pinning sites and the resulting j c {B) consists of all the different interactions. It 
is relatively simple to calculate the pinning force of a single type of pinning site at low 
field^, but even increasing the field, when vortex-vortex interactions come into play, makes 
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the calculations in closed form impossible. Statistical approach has been used e.g. in refs. 



19-22 



and these models describe usually the shape of the experimental curves very well, 
but the understanding of e.g. the vortex paths inside the superconductors is limited. 

The Ginzburg-Landau equations, although originally meant to be used close to the phase 
transition, describe superconductors well at lower temperatures too. Solving the minimum 
energy configuration for a certain set of parameters gives us the spatial variation of the order 
parameter ip, which can be used to calculate the fraction of the pinned vortices as a function 
of e.g. the magnetic field or the density, shape and size of the pinning sites. The Ginzburg- 
Landau equations have been solved for small particles^ and with one or several pinning 
sites using 5T C pinning^ and by restricting the order parameter—. However, Crabtree et a/.- 
present a large scale model with flux pinning that is similar to this work. Unfortunately, 
they do it for a relatively low k — 4 (low for modelling YBCO) and concentrate on the 
dynamics of the vortex trajectories instead of the critical current that is the focus of this 
work. We present a method for doing large scale computation of pinning in superconductors 
that are close to realistic size and calculate the field depedence of the critical current. In 
this paper we consider the case where magnetic field is parallel to the columnar defects. The 
results are directly compared to data from thin film YBCO samples. 



II. COMPUTATIONAL MODEL 

A. Solving the Ginzburg-Landau equations 

We chose to solve the static Ginzburg-Landau equations by finding a (local) minimum 
of the associated energy functional. For computational purposes, we write the energy in a 
dimensionless form. The only dimensional value is the overall energy scale, which does not 
affect the solutions and the dimensionless energy is 

E= /"d 3 x(|||(V + zi)^|| 2 + |||Vxl|| 2 

J (1) 

+X(h/f-i) 2 ), 

where k = y^5/ (2fi h 2/ y 2 q 2 ) is the dimensionless Ginzburg-Landau parameter and 7 = 
l/(4m e ) and q = 2e, the penetration depth is absorbed in the overall energy scale and 
therefore A = 1. Naturally, the coherence length is now simply £ = 1/k. Equation ([T|) 
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is then discretised as described in ref. |26[ In short, this discretization preserves the gauge 
invariance of the system and allows us to solve the equations without choosing and enforcing 
a gauge and the associated problems. 

The solutions are then found using the TAO^Z and PETSo^ - — massively parallel nu- 
merical libraries. We account for the impurities in the physical system by using a bound 
constrained variant of the limited memory quasi-Newton algorithm (also called a variable 
metric algorithm) with BFGS^~— formula for Hessian approximations. This is provided by 
the TAO library. The impurities are then modelled by setting appropriate constraints for 
the Cooper pair density at the impurity locations. The gradients required by the algorithm 
are computed from the discretized energy in a straightforward manner. 

The correctness and accuracy of the program was tested by solving the equations for a 
single Abrikosov vortex. At a 102 2 lattice and lattice constant of 0.1 ~ 0.07£, the deviation 
from correct total energy is less than 0.02 % with most of the deviation resulting from 
the small box: increasing the number of lattice points increases the accuracy. A smaller 
lattice constant does not increase the accuracy as much. Note that due to the rescaling, 
the lattice constant is in units of penetration depth A. To be able to resolve the core of the 
vortex the lattice constant has to be smaller than £. This means that a lattice constant h 
of approximately 0.1£ gives very accurate results. Computer memory is a limiting factor in 
the computations. In practice, h ~ 0.3£ yields accurate enough results and requires ~ 97 % 
less computer memory. 

At the y- and ^-boundaries of the calculation ijj is set to zero and the vector potential is 
kept fixed to simulate an external magnetic field. The x-boundary is periodic. Lattice sizes 
used in simulations were typically 500 x 520 x 50 of which 10 pixels near the boundaries were 
used as vacuum (ip = 0, A is free). The vacuum between the sample and the calculation 
boundary allows for the magnetic field to bend around the sample which makes the external 
field of the simulation comparable to the actual field of the measurement of a thin film in a 
magnetometer. 

As an example of a simulation result fig. [1] shows the absolute value of i/j (fig. (IK) and the 
phase (fig. [Tb) with B = 3 T in a sample with dislocations. Using also the vector potential 
A, the magnetic field B z = (V x A) z (fig. [Tb) and the current in the periodic direction 
j x = (V x B) x /fi (fig. [TH) were calculated. Note that also the shielding currents are visible. 
The calculation always yielded zero current in the vacuum, thus confirming the validity of 



the vacuum spacer layer. 



B. Introducing pinning to the model and calculating j c 

Flux pinning was modelled by locally restricting the maximum possible value of the real 
and imaginary parts of ip which then represents a pinning site. The chosen maximum value 
was 0.1 which corresponds to limiting the maximum value of \ip\ between 0.1 and a/2/10. 
The limit was not set to zero because that would have made the analysis of the result 
considerably more difficult. With this limit vortices can be defined to be the region of the 
calculation lattice where < 0.1. 

The pinning sites had the same physical size and shape as is deduced from transmission 
electron microscopy (TEM) data: The BZO-nanorods were modelled as randomly distributed 
sample penetrating rods with a diameter of 5 nm^i^. The dislocations were c-axis aligned 
sample penetrating rods with a diameter of 0.3 nm^, which were also randomly distributed 
in the sample. 

The critical current was determined as the fraction of the vortex length trapped in the 
pinning site to the total length as in ref. O. Thus, we do not get the absolute values, but 
the field dependence of the j c instead, which can be directly compared to the experimental 
data. The total length of vortices was calculated by following each vortex. If a point along 
the vortex is closer than 1.5£ to a pinning site it was counted to the pinned section of the 
vortex. Finally, j c is proportional to the fraction of the pinned sections to the total vortex 
length. 

Modelling different pinning site types at the same time requires calculations over different 
length scales which is memory consuming. Dislocations are small in size and sparsely dis- 
tributed which means large sample sizes are needed to contain more than a few dislocations. 
As a solution to this memory issue the coherence length and the pinscape were scaled up so 
that k = 10. Thus, we can use larger lattice constant allowing sparse dislocation densities. 
The change of k can be done without major changes to the physics^, since even at k — 10, 
we are at the limit of high-K and the magnetic field variation inside the superconductor is 
small. This was also verified by calculating two simulations with the same pinscape and 
resolution in units of £ but with different n. Figure [2] shows the results of such calculations. 
It is easily seen that changing k does not have an appreciable effect on the results. 
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III. RESULTS AND DISCUSSION 



The field dependence of j c with BZO-nanorods as a pinscape was simulated for a 230 nm 
wide sample with thickness of 15 nm that is thick enough for the order parameter to reach 
1 inside the sample and for all the relevant physics since here B is always perpendicular to 
the sample. The third direction was periodic with a period of 240 nm. There was also a 5 
nm thick layer of vacuum around the sample. Both the rods and the magnetic field were 
parallel to c-axis. The magnetic field was varied from 4 T to T with small steps. The 
result of the calculation with the previous field value was used as the initial condition for 
the calculation with the next field value. The density of the BZO-nanorods was set so that 
the matching fields were = 1 T (30 rods) and B^ = 2 T (60 rods). 

The j c with dislocations as pinning sites was calculated with a similar calculation grid but 
its length was scaled by a factor of ten, due to the change in k, to achieve the required low 
density of dislocations. The fine tuning of the dislocation density to match the experimental 
data was done after the simulation by scaling the length of the calculation lattice unit cell. 
The scaling does not alter the results because the pinning strength is determined by the 
ratio r r /£, and vortex- vortex interactions are affected by the ratio of the average pinning 
site separation to the penetration depth A which is also dimensionless. What the scaling does 
affect are the size of the sample and the value of the magnetic field. Thus, the sample size 
for the dislocation simulations was 530 nm (periodic) x 530 nm x 34 nm with dislocation 
densities corresponding to matching fields of 90 mT (12 rods), 180 mT (25 rods) and 360 
mT (50 rods). 

Figure [3] shows examples of simulation results at fields B = 3 T and B = 1.5 T with 
BZO-rods and dislocations as pinning sites. It can be seen that the strong pinning force of 
BZO-nanorods strongly disturbs the vortex lattice while pinning by dislocations results in a 
more regular vortex lattice. The modelled pinning sites can be seen in these images as larger 
circles (nanorods) and small dots (dislocations). It is easy to determine whether a pinning 
site is occupied when these images are combined with the information about the phase of i/j 
and the pinscape coordinates. 

The magnetic field scan from high field to low field was calculated using several random 
pinscapes with the same number of pinning sites. The results were taken as the averages of 
the critical currents derived from these simulations and the errors as the standard deviations 
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of the simulation sets. In figure H] we have overlayed the obtained j c {B) for nanorods with 
2.9 wt-% and 9 wt-% BZO-doped YBCO film data±± and the j c (B) for dislocations with 
data from an undoped YBCO film. The j c values have been scaled to 1 at zero field. 
The agreement with the experimental data is excellent and shows that the nanorods are 
so dominant in pinning that other pinning types can be neglected when the rods and the 
magnetic field are parallel. The a values for dislocation simulations range from 0.5 to 0.8 
while for BZO-nanorods simulations gave a values of 0.26 and 0.41 for rod densities of 1 T 
and 2 T, respectively, which are quite typical v&luesM-^^— . The a value of the measured 
pure YBCO thin film shown in fig. Sb is 0.59. Comparing the dislocation data to the 
experiments is more challenging because experimentally dislocation density is not so easy to 
adjust with the growth conditions while the density of the nanorods can be adjusted with 
the doping level. Further difficulties arise from the fact that dislocations do not dominate 
pinning so clearly as the large nanorods do. Thus, including other types of pinning sites (e.g. 
twins, oxygen vacancies) to the model is needed for more in depth analysis of the undoped 
films. 

The nanorod data in fig. [4^ shows a rounded shape of j c {B) on log- log scale which is 
in contrast to the sharp bend between the different field regimes of the dislocation data in 
fig. Hb. Fitting B~ a to the rounded curves is difficult because the shape of the curve is not 
really correct unlike with dislocations. In this work a was determined from the part close 
to B* if the curve is rounded. 

It would seem obvious to attribute the rounded shape of the BZO-nanorod data to the 
more dominant vortex-vortex interactions caused by the higher density of the defects. But 
the dislocations differ from nanorods also by size not only by density. Thus, we calculated 
j c for a high density (B^ = 1.5 T) pinscape with fixed defect locations but with variable rod 
radius r r which is shown in fig. From this it is clear that the typical rounded shape of 
the j c -curve of the BZO-nanorods requires not only a high density but also a large rod size. 
Having a high density of dislocations will not make the sample as good as one with a high 
density of BZO-nanorods. 

The results in fig. were further analyzed by fitting j c (B) = cB~ a to the above B* 
portion of the datapoints. The obtained a values are shown in fig. The values decrease 
with increasing rod diameter from 0.6 to 0.2 which is within the range of measured values 
for BZO and BSO nanorods schematically shown as ellipses in fig. [6j Naturally, the decrease 
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of a comes from the deeper pinning potential of a larger rod which makes it possible for a 
vortex to sit in the potential well up to higher magnetic fields. The ellipses fit the simulation 
results even better if we take into account that in real samples nanorods are surrounded by 
additional defects which makes the rod effectively larger than the actual rod size measured 
with TEM and would thus move the ellipses to the right. 

Furthermore, the accommodation field was calculated from the fit as B* = %/c i.e. the 
point where cB~ a intersects the line j c {0) = 1. The accommodation fields relative to the 
matching field are shown in fig. [7J The accommodation fields level off to value B* / = 0.7 
starting from the rod size r r /£ = 2 which is the size where several vortices get pinned at a 
single rod. The B* /B^ = 0.7 is also seen in experiments^*^. Assuming strong pinning by 
linear defects in low magnetic field the pinning potential per unit length e r of the vortex is 
related to the accommodation field^: 

B e r . „ 

where e = irk 2 /(fioq 2 \ 2 ) is the characteristic energy of the vortex per unit length. An 
analytical expression for the depth of the pinning potential of a cylindrical cavity has been 
derived in ref. I^l as an upper limit within the London approximation: 



Using eqs. [2] and [3] the accommodation field can be roughly related to the pinning potential 
of the nanorods. But as can be seen from fig. [7J this relation between the accommodation 
field and the pinning potential seems to break down at large rod sizes where several vortices 
are pinned per rod. 

The pinning potential of a nanorod was also determined by simulating a system consisting 
of a single vortex and a single nanorod. The depth of the potential was taken as the difference 
in the total energy between the state where the vortex is far away from the pinning site and 
where the vortex is pinned in the nanorod divided by the length of the rod. This was done 
for several rodsizes and the result is shown in fig. [7J where eq. [3] and the accommodation 
fields are also shown. Considering that there is no scaling or fitting involved the calculated 
potential depths are in general agreement with eq. [3J 

Examples of more than one vortex pinned to a nanorod are shown in fig. [HI Up to four 
vortices get pinned to a large nanorod at high magnetic field. The vortices sit symmetrically 




(3) 
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at opposite sides of the rod. At rod size r r /£ = 2 the vortices are forced to be very close 
to each other but still there are two vortices in two of the rods. At rod size r r /£ = 1 there 
is only one vortex per rod. Ideally, there should be one vortex at each pinning site at the 
matching field. But the pinning sites are randomly distributed which makes them unevenly 
spaced. Thus, at small rod sizes the free vortices and the unoccupied pinning sites cancel out 
each other in the large scale giving the average of one vortex per pinning site. At large rod 
sizes this cancelling out happens by pinning several vortices at suitable sites while leaving 
some sites empty. 

It is obvious from fig. |H] that the vortex lattice is more regular with the small nanorods. 
For a closer look into the short and long range order the radial distribution functions (fig. 
[9]) of the vortex positions (N ~ 250) were calculated from the simulation results at 7.5 T. 
The vertical lines show the positions of the maxima for an ideal 2D-triangular lattice with 
the unit length of a = (2$o/("\/3-B)) 1 ^ 2 where $o is the flux quantum and B = 7.5 T. There 
are 4 broad peaks visible at the rod size 0.3 which roughly overlap the positions of the ideal 
case. The 4th peak corresponding the 6th and 7th nearest neighbour (NN) is already very 
weak. There is no long range order past the 7th NN at any rod size. At the rod size r r /£ = 3 
and larger only the first NN peak is visible. Thus, the triangular vortex lattice continues 
over the small rods with some disturbance while the large rods completely break down the 
long range order. Since the area of the first peak is constant the coordination number is 
the same for all rod sizes. The pinning of several vortices to the same pinning site has been 
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where the irregularity of the vortex lattice near the 



experimentally observed e.g. in ref. 
pinning sites is clearly visible too. 

From the experimental point of view the results presented here further emphasize that one 
should not only focus on optimizing the density of the pinning sites but the size of the defects 
needs to be carefully considered too. Naturally, if the magnetic field in the application is not 
homogenous one needs to consider the splay of the nanorods too. In large scale production, 
where in situ deposition is almost the only option, the size of the defect is fixed by the 
chemical properties of the dopant^. Thus, one has to try different dopant materials to find 
those that produce large enough columnar defects. With large defects the optimal density 
is a compromise between the superconducting volume and pinning. However, at large defect 
sizes, where several vortices get pinned to each defect, it is clear that any density (in terms 
of matching field) of the defects above the operating field of the application is not optimal 



since some of the defects will be left empty. On the other hand, if there are two vortices per 
rod, a density of half the operating field is too low because some of the rods will only pin 
one vortex due to e.g. the variation in the rod positions. With naive arguments one can say 
that having one vortex per rod is better than two since with two vortices the rod diameter 
has to be roughly twice as large which increases the volume of the rods by a factor of 4. 
But in reality there is a lot of different strains and relaxations involved when the distance 
between the rods is changed which could compensate for the loss of the superconducting 
volume by increasing T c and j c (0). 

IV. CONCLUSIONS 

In this paper we have shown that it is possible to model flux pinning numerically with 
the Ginzburg-Landau equations. The model fully includes vortex-vortex interactions which 
is very important for modelling flux pinning. Doing this allows us to see the vortex paths 
inside the sample and to calculate the field dependence of j c for a pre-determined pinning 
site configurations. The results are in excellent agreement with the experiments. 

The main result was that the reason behind the characteristic round shape of the j c and 
the lowered a value of the BZO-doped YBCO films is not only the high density but also the 
large size of the pinning sites. The density of the pinning sites acts through the vortex- vortex 
interactions while the rod size (i) changes the depth of the pinning potential, (ii) changes 
the number of the pinned vortices per rod, (iii) and supresses the long range order of the 
vortex lattice at larger sizes. 

The current model can also be used to simulate the angular depedency of j c which will 
be our next objective. Adding time and current to the model by using time-dependent 
Ginzburg-Landau will allow transport measurement simulations. 
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FIGURES 
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FIG. 1. a) The absolute value of the calculated order parameter, shown on a plane sliced 
perpendicular to c-axis. b) The phase of the order parameter in the same calculation, c) The 
magnetic field density calculated from the vector potential in the same simulation and d) the 
current flowing in the periodic direction (up and down) . Note that the scale is cut off short so that 
also the shielding current is visible. 
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FIG. 2. The same simulation calculated with k = 10 and k = 100. The change in k does not 
change the results. 
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FIG. 3. The absolute value of the calculated order parameter, \ip\, shown on a plane sliced 
perpendicular to c-axis for nanorods (on the left) and dislocations (on the right) as pinning sites. 
The applied fields are 3 T (top) and 1.5 T (bottom). The scale of the images is the same as in fig. 
Eh, 
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FIG. 4. Calculated j c (empty symbols) compared to experimental data (linespoints) when the 
simulation and sample contain a) nanorods and b) dislocations. The simulation point is an average 
over several different pinscapes and the errorbars are the standard deviations of those. 
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FIG. 5. Calculated j c with the same rod locations but different rod radii r r . 
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FIG. 6. The a values from fits to curves in fig. [5j The error bars are the standard deviations of 
the fits. The range of experimental values measured for pur o 36l37 ; BZO doped^— and BSO*^i 
doped YBCO thin films is indicated with ellipses. 
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FIG. 7. The depth of the pinning potential 4e r /eo as a function of the rod radius r r as calculated 
from the simulations (blue squares) and as given by eq. [3] (solid curve). The red circles show the 
ratio of the accommodation field B* to the matching field determined from the data shown in 

fig. El 




FIG. 8. A closer look at 3 of the 55 nanorods in the simulations with varying rod size. The 
absolute value of the order parameter is shown with the color scale from 1 (blue) to (red). The 
size of the images is 230 x 160 simulation grid points (130 nm x 90 nm). Simulation results in 
high magnetic field are on the left and at matching field (1.5 T) are on the right. The size of the 
nanorod decreases from top to bottom. 



17 




10 20 30 40 50 60 70 
Distance (nm) 



FIG. 9. The radial distribution functions of the vortex core positions for the different nanorod 
sizes. The vertical lines mark the positions of the nearest neighbours for an ideal 2D-triangular 
vortex lattice at 7.5 T. 
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